Polyfem Examples¶
This is a jupyter notebook. The "real" notebook can be found here.
Some imports for plotting
import plotly.offline as plotly
import plotly.graph_objs as go
#Necessary for the notebook
plotly.init_notebook_mode(connected=True)
algebra
import numpy as np
stuff for the animation
import ipywidgets as widgets
from ipywidgets import interact
and finallypolyfempy
import polyfempy as pf
Plotting utility¶
This code is not particularly interesting.
It converts a tet-mesh (4 indices) into a face-based tri mesh and uses plotly Mesh3d to plot it.
def plot(vertices, connectivity, function):
x = vertices[:,0]
y = vertices[:,1]
if vertices.shape[1] == 3:
z = vertices[:,2]
else:
z = np.zeros(x.shape, dtype=x.dtype)
if connectivity.shape[1] == 3:
f = connectivity
else:
#Convert a tet-mesh into face based triangles.
f = np.ndarray([len(connectivity)*4, 3], dtype=np.int64)
for i in range(len(connectivity)):
f[i*4+0] = np.array([connectivity[i][1], connectivity[i][0], connectivity[i][2]])
f[i*4+1] = np.array([connectivity[i][0], connectivity[i][1], connectivity[i][3]])
f[i*4+2] = np.array([connectivity[i][1], connectivity[i][2], connectivity[i][3]])
f[i*4+3] = np.array([connectivity[i][2], connectivity[i][0], connectivity[i][3]])
mesh = go.Mesh3d(x=x, y=y, z=z,
i=f[:,0], j=f[:,1], k=f[:,2],
intensity=function, flatshading=function is not None,
colorscale='Viridis',
contour=dict(show=True, color='#fff'),
hoverinfo="all")
layout = go.Layout(scene=go.layout.Scene(aspectmode='data'))
fig = go.Figure(data=[mesh], layout=layout)
plotly.iplot(fig)
Set the mesh path
mesh_path = "plane_hole.obj"
create settings
settings = pf.Settings()
pick linear $P_1$ elements. If the mesh would be a quad it would be $Q_1$
settings.discr_order = 1
normalize the mesh to be in [0,1]^2
settings.normalize_mesh = True
and choose Young's modulus and poisson ratio
settings.set_material_params("E", 210000)
settings.set_material_params("nu", 0.3)
We are use a linear material model
settings.tensor_formulation = pf.TensorFormulations.LinearElasticity
Next we setup the problem
problem = pf.GenericTensor()
sideset 1 has zero displacement in $x$
problem.add_dirichlet_value(1, [0, 0], [True, False])
sideset 4 has zero displacement in $y$
problem.add_dirichlet_value(4, [0, 0], [False, True])
sideset 3 has a force (Neumann) of [100, 0] applied
problem.add_neumann_value(3, [100, 0])
fianally set the problem
settings.set_problem(problem)
Solve!
solver = pf.Solver()
solver.settings(str(settings))
solver.load_mesh(mesh_path)
solver.solve()
Get the solution
[pts, tets, disp] = solver.get_sampled_solution()
diplace the mesh
vertices = pts + disp
and get the stresses
mises, _ = solver.get_sampled_mises_avg()
finally plot with the above code
plot(vertices, tets, mises)
Note that we used get_sampled_mises_avg to get the Von Mises stresses because they depend on the Jacobian of the displacement which, becase we use $P_1$ elements, is piece-wise constant. To avoid that effect in get_sampled_mises_avg the mises are averaged per vertex weighted by the area of the triangles. If you want the "real" mises just call
mises = solver.get_sampled_mises()
plot(vertices, tets, mises)